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Abstract 

We present new iterative algorithms for solving a square linear system Ax = b in dimension n by 
employing the Triangle Algorithm [IT] . a fully polynomial-time approximation scheme for testing if the 
convex hull of a finite set of points in a Euclidean space contains a given point. By converting Ax = b 
into a convex hull problem and solving via the Triangle Algorithm, together with a sensitivity theorem, 
we compute in 0(n 2 e~ 2 ) arithmetic operations an approximate solution satisfying \\Ax e — b\\ < ep, 
where p = max{||ai||, . . . , ||a n ||, ||6||}, and a» is the i-th column of A. In another approach we apply the 
Triangle Algorithm incrementally, solving a sequence of convex hull problems while repeatedly employing 
a distance duality. The simplicity and theoretical complexity bounds of the proposed algorithms, requiring 
no structural restrictions on the matrix A, suggest their potential practicality, offering alternatives to 
the existing exact and iterative methods, especially for large scale linear systems. The assessment of 
computational performance however is the subject of future experimentations. 

Keywords: Linear System of Equations, Iterative Methods, Convex Hull, Linear Programming, Duality, 
Approximation Algorithms. 

1 Introduction 

The significance of methods for solving a linear system of equations becomes evident with high school 
mathematics and science classes. Solving linear system of equations is undoubtedly one of the most practical 
problems in numerous aspects of scientific computing. Gaussian elimination is the most familiar method for 
solving a linear system of equations, discussed in numerous books, e.g. Atkinson [T], Bini and Pan [3J, Golub 
and van Loan |13j , and Strang |24] . The method itself is an important motivation behind the study of linear 
algebra. Iterative methods for solving linear systems offer very important alternatives to direct methods and 
find applications in problems that require the solution of very large or sparse linear systems. For example, 
problems from discretized partial differential equations lead to large sparse systems of equations making the 
use of direct methods impractical. 

Iterative methods generate a sequence of approximate solutions. They begin with an initial approxi- 
mation and successively improve it until a desired approximation is reached satisfying a certain measure 
of error. Iterative methods include such classical methods as the Jacobi, the Gauss-Seidel, the successive 
over- relaxation (SOR), and the symmetric successive over-relaxation method which applies to the case when 
the coefficient matrix is symmetric. When the coefficient matrix is symmetric and positive definite the con- 
jugate gradient method (CG) also becomes applicable. Convergence rate of iterative methods can often be 
substantially accelerated by preconditioning. Some major references in the vast subject of iterative methods 
include, Barrett et al. [2 , Golub and van Loan [13], Greenbaum [15], Saad [23], van der Vorst [25], Varga 
[27] . and Young [28] . 

To guarantee convergence, iterative methods often require the coefficient matrix to satisfy certain condi- 
tions, or decompositions, necessary to carry out the iterations. Some of the theoretical analysis, especially 
in the earlier works on iterative methods, is only concerned with convergence or conditions on the coefficient 
matrix or its decompositions that guarantee convergence. However, in some cases theoretical complexity 



analysis is provided, see e.g. Reif [32] who considers the complexity of iterative methods for sparse diago- 
nally dominant matrices. 

The major computational effort in each iteration of the iterative methods involves matrix-vector multi- 
plication. This makes iterative methods very attractive for solving large systems and also for parallelization. 
There is also a vast literature on parallelization of iterative methods for large systems, see Demmcl[6 , 
Dongarra et al. [7] , Duff and van der Vorst [8] , van der Vorst [35] , van der Vorst and Chan [3j5] . 

In this article we offer a very different iterative method for solving an n x n linear system Ax = b. 
It is inspired by a new geometric algorithm for a convex hull problem in |17j . To illustrate the approach 
to be presented, we begin with a very simple example. Consider a triangle with vertices v\ = (wn,Ui2), 
«2 = («2ij V22)) ^3 = (^31,^32), and a point p = {p\,P2) in their convex hull (Figure [1]). 




Figure 1: A point p in the convex hull of vertices of a triangle. 



In particular, there exists nonnegative scalars ati, 0*2,0:3 such that 

aivi + a 2 v 2 + CX3V3 = P, 011 + a 2 + a 3 = 1. (1) 

When the vertices are not collinear, the scalars can be determined by solving a linear system having an 
invertible coefficient matrix: 

(V\\ V21 V3i\ I Ctl\ I Pl\ 

(2) 

However, since p is known to lie in the convex hull of the vertices, we can solve this system using a new 
geometric approach. In this article we will first review our iterative algorithm for solving the general case of 
this convex hull problem. Then, we will show how to convert the general linear system Ax = b into a convex 
hull problem so that the convex hull algorithm would be applicable. 

In general the convex hull of a set S in a Euclidean space is the intersection of all convex sets that contain 
it. In particular, if S — {vi, . . . , v n } C K m , its convex hull written as conv(S) is given by 




i(S) = -j am, } aj = l, ai > 0, Vi = l,...,n> 

^ i=l i=l > 



(3) 



In this article we describe new iterative methods for solving a linear system of equations utilizing a new 
algorithm, called the Triangle Algorithm [17] . designed to solve the following convex hull decision problem: 

Given a set of points S = {vi, ■ ■ ■ , v n } C M. m and a distinguished point p G R m , test if p lies in conv(S). 
A practical approximate version of the convex hull decision problem is the following: 

Given e £ (0, 1), either compute a point p e € conv(S) such that 

\\p - Pe\\ < e\\p - Vi\\, for some i = 1, n; (4) 

or prove that p ^ conv(S). The complexity of the Triangle Algorithm in computing an approximate solution 
as in (j4]), assuming such solution exists, is 0(mne~ 2 ) arithmetic operations. This complexity is attractive in 
contrast with the complexity of algorithms that solve the convex hull problem exactly. 

The convex hull decision problem is a basic and fundamental problem in computational geometry as well 
as in linear programming (LP). More general convex hull problems are described in Goodman and O'Rourke 
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P3] ■ On the one hand the convex hull decision problem is a very special case of LP feasibility. On the other 
hand, it is well known that through LP duality theory the general LP problem may be cast as a single LP 
feasibility problem, see e.g. Chvatal [4 . This problem in turn is related to the convex hull decision problem. 
In fact it can be justified that the two most famous polynomial-time LP algorithms, the ellipsoid algorithm 
of Khachiyan [115] and the projective algorithm of Karmarkar [3T] , are in fact explicitly or implicitly designed 
to solve a case of the convex hull decision problem where p = 0, see |16j . In |20) another polynomial-time 
algorithm is given for this special homogeneous case based on diagonal matrix scaling. For integer inputs all 
known polynomial-time LP algorithms would exhibit theoretical complexity that is polynomial in m, n, and 
the size of encoding of the input data, often denoted by L, see e.g. [TH]. The number L can generally be 
taken to be dependent on the logarithm of ran and of the logarithm of the absolute value of the largest entry 
in the input data. No strongly polynomial-time algorithm is known for LP, i.e. an algorithm that would in 
particular solve the convex hull decision problem in time complexity polynomial in ra and n. 

The convex hull decision problem can also be formulated as the minimization of a convex quadratic 
function over a simplex. This particular convex program has found applications in statistics, approximation 
theory, and machine learning, see e.g Clarkson [3] and Zhang [29] who consider the analysis of a greedy 
algorithm for the more general problem of minimizing certain convex functions over a simplex (equivalently, 
maximizing concave functions over a simplex) . The oldest version of such greedy algorithms is Frank- Wolfe 
algorithm, [9]. Special cases of the problem include support vector machines (SVM), and approximating 
functions as convex combinations of other functions, see e.g. Clarkson [5]. The problem of computing the 
closest point of the convex hull of a set of points to the origin, known as polytope distance is the case of the 
convex hull decision problem where p is the origin. In some applications the polytope distance refers to the 
distance between two convex hulls. Gilbert's algorithm [TT], [T3] for the polytope distance problem is one of 
the earliest known algorithms for the problem. Gartner and Jaggi [10] show Gilbert's algorithm coincides 
with Frank- Wolfe algorithm when applied to the minimization of a convex quadratic function over a simplex. 
In this case the algorithm is known as sparse greedy approximation. 

From the description of Gilbert's algorithm in [10] it does not follow that Gilbert's algorithm and the 
Triangle Algorithm are identical. However, there are similarities in theoretical performance of the two algo- 
rithms and in [17] we give such theoretical comparisons between the Triangle Algorithm for the convex hull 
decision problem and the sparse greedy approximation. Indeed we believe that the simplicity of the Triangle 
Algorithm and a new duality theorem that inspires the algorithm, as well as its theoretical performance all 
make it distinct from other algorithms for the convex hull decision problem. Furthermore, these features of 
the Triangle Algorithm may encourage and inspire new applications of the algorithm. The present article is 
testimonial to this claim where we will use the algorithm to solve a linear system. 

In what follows we present two iterative algorithms for solving an invertible matrix equation Ax = 
b in dimension n via the Triangle Algorithm. In the first algorithm we convert Ax — b into a convex 
hull problem. Then given e > 0, we apply the Triangle Algorithm and a sensitivity theorem from [17] 
to compute in 0(n 2 e~ 2 ) arithmetic operations an approximate solution satisfying \\Ax e — b\\ < ep, where 
p = max{||&||, ||ai||, . . . , ||a„||}, and is the i-th column of A. In the second algorithm we apply the Triangle 
Algorithm incrementally, solving a sequence of convex hull problems while using a distance duality from [17] . 
The first algorithm requires an a priori estimate of the least coordinate of the solution while the second 
algorithm builds up this estimate iteratively. These offer alternatives to exact methods and to iterative 
methods such as the Jacobi, Gauss-Seidel, and successive over-relaxation methods, requiring none of their 
structural restrictions. The extreme simplicity and theoretical complexity bounds suggest the potential 
practicality of the methods, especially for large scale or sparse systems. The assessment of computational 
performance is the subject of future experimentations. In 18 we considered an application of the Triangle 
Algorithm for solving a very special case of a matrix equation where the vector b is unknown but some 
attributes such as its Euclidean norm is available. 

This article is organized as follows. In Section [5] we describe the Triangle Algorithm from [T7J , a duality 
that it makes use of, the distance duality, as well as some geometric properties of this duality. In Section 
[3] we consider solving a linear system, Ax — b. In 13.11 we consider the case where x = A~ 1 b is known to 
be nonnegative and show how the Triangle Algorithm together with a sensitivity theorem from [17j solves 
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a linear system and state its complexity for computing an approximate solution. In 13.21 we consider the 
general case of Ax = b where there is no information on the solution and show how this can be converted 
into the first case by computing an a priori lower bound t* on the minimum component of the solution 
vector x — A~ 1 b. Then we state a complexity result. In 13.31 we offer an incremental version of the Triangle 
Algorithm for solving Ax — b not requiring an a priori lower bound on t*. We conclude with some final 
remarks. 

2 Review of The Triangle Algorithm 

Throughout we let || • || denote the Euclidean norm. In |17j we proved the following characterization theorem 
for the convex hull decision problem. 

Theorem 1. (Distance Duality [T7] ) Let S = {vi, . . . , v n } C R"\ p £ R m . Then we have 

(i) : p £ conv(S) if and only if given any p' £ conv(S), there exists Vj such that \\p' — Vj\\ > \\p — Vj\\. 

(ii) : p £ conv(S) if and only if there exists p' £ conv(S) such that \\p' — V{\\ < \\p — Vi\\, V i. □ 

Definition 1. We say p' £ conv(S) is a witness if it satisfies \\p' — fill < \\p — Vi\\, for all i = 1, . . ., n. 

Each witness certifies the infeasibility of p in conv(S). The next theorem shows that each witness actually 
induces a separating hyperplane. The set W p of all such witnesses is the intersection of conv(S) and open 
balls Bi of radius \\p — Vi\\ centered at Uj, i = 1, . . . , n and forms a convex open set in the relative interior of 
conv(S). 

Theorem 2. (Characterization of Witness Set [17J) p' £ W p if and only if the orthogonal bisector 
hyperplane of the line segment pp' separates p from conv(S). □ 



Figure 2: Examples of empty W p (p £ conv(S)) and nonempty W p (p & 1 conv(S)), gray area. 

The following is a simple but useful consequence of Theorem [2] 
Corollary 1. Suppose p ^ conv(S). Let 




A = min{||p — x\ 



x £ conv(S)}. 



(5) 



Then any p' £ W p gives an estimate of A to within a factor of two. More precisely, 




(6) 
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Definition 2. Given e <E (0, 1) we say p' € conv(S) is an e-approximate solution if it satisfies 



\\p'-p\\<eR, (7) 

where 

R = max{\\p-v 1 \\,...,\\p-v n \\}. (8) 

Using the characterization theorem (Theorem [lj , in [17] we described a simple algorithm, called the 
Triangle Algorithm. Given a desired tolerance e € (0, 1), and a current iterate p' £ conv(S), in each iteration 
the Triangle Algorithm searches for a triangle App'vj where Vj € S 1 satisfies — > \\p — Vj\\. Given that 
such triangle exists, the algorithm uses Vj as a pivot point to "pull" the current iterate p' closer to p to get 
a new iterate p" € conv(S). If no such a triangle exists, then by Theorem[Tl p' is a witness certifying that p 
is not in conv(S). The Triangle Algorithm consists of iterating two steps: 



Triangle Algorithm (S — {vi, . . . , v n }, p) 

• Step 1. Given an iterate p' € conv(S) \ {p}, check if there exists a pivot point Vj € S, i.e. 
\\p' — v j\\ > \\p — Vj\\. If no such Vj exists, then p' is a witness, stop. 

• Step 2. Otherwise, compute the step-size 

_ {p-p') T {v 3 -p') 



Let the iterate be defined as 



(1 — a)p' + avj, if a £ [0, 1]; 



pll= ^ „„„ (10) 



Replace p' with p", go to Step 1. 



I vj , otherwise 



By an easy calculation that shift p' to the origin, it follows that the point p" in Step 2 is the closest point 
to p on the line p'vj , see Figure [3] 




Figure 3: Depiction of gaps S = d(p',p), 5' = d(p",p). 



Since p" is a convex combination of p' and Wj, it will remain in conv(S). The algorithm replaces p' with 
p" and repeats the iterative step. Note that a pivot point Vj may or may not be a vertex of conv(S). In fact 
p" can be explicitly written as a convex combination of v^s, 

n 

p" = ^2fiiVi, /3j = (1 - a)aj + a, (3, t = (1 - a)a i7 Vi ^ j. (11) 

i=l 

The following straightforward result implies that testing for a pivot point or a witness we do not need to 
compute square-roots. 
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Proposition 1. A point Vj is a pivot point if and only if 



(P-P') T VJ>1(\\P\\ 2 -\\PT)- (12) 



Equivalently, p' is a witness if and only if 



(p-p') T « i <i(||p|| 2 -||p'|| 2 ), Vi = !,...,». (13) 



Remark 1. If p = 0, then Vj is a pivot point if and only if 

P ,T v 3 < l\\p'\\ 2 - (14) 
In particular, if p' T Vj < 0, then Vj is a pivot point. Equivalently, p 1 is a witness if and only if 

P' T ^ > \\\p'\\\ Vi = l,...,n. (15) 

The following theorem in [T7] estimates the complexity of the Triangle Algorithm. 

Theorem 3. The Triangle Algorithm correctly solves the convex hull problem as follows: 

(i) Suppose p G conv(S). Given e > 0, the number of iterations K e to compute a point p e in conv(S) so 
that 

\\Pe~p\\ < 4p~ v j\\> (16) 

for some Vj G S satisfies 
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K t <-=0{e-'). (17) 
e £ 

(ii) Suppose p conv(S). If A — min |||a; — p\\ : x G conv(S)}, the number of iterations to compute 
a witness, a point p a in conv(S) so that d{p/^ 1 Vi) < d(p,Vi) for all Vi G S, satisfies 



K A < 



8R\ (25 



□ (18) 



3 Solving A Linear System Via The Triangle Algorithm 

Consider solving Ax = b with A an invertible n x n real matrix. Let A — [ai, G2, • • • , a n ], where a, G K™. 
Definition 3. We say xq is an e-approximate solution of Ax = b if 

\\Ax - b\\ < ep, p = max{|| ai ||,...,||a„||,||6||}. (19) 

3.1 Solving A Linear System With Nonnegative Solution 

First, we assume that x = A~ 1 b > and show how to solve this as a convex hull problem. Next, we solve the 
general case relaxing this condition. Since A is invertible, Ax = has only the trivial solution. In particular, 

gconv({ai,...,a n }). (20) 

In [17] we described an application of the Triangle Algorithm in computing an approximate feasible 
solution of the feasibility problem in linear programming. Here we use that approach to solve Ax = b via 
the Triangle Algorithm to give for any e G (0,1), an e-approximate solution. Clearly, by adjusting e we 
can compute an approximate solution with absolute error of e, however the approximate solution within a 
relative error as defined in Definition [3] is a more sensible measure of approximation. The following auxiliary 
result is easy to prove. 
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Proposition 2. Assume x = A 1 b>0. Then 6 conv({ai, . . . , a n , — b}). □ 

It follows that solving Ax — b approximately is equivalent to finding an appropriate approximation to 
in the set conv({ai, . . . , a n , —b}). 

The following sensitivity theorem establishes the needed accuracy to which an approximate solution in 
conv({ai, . . . , a n , —b}) should be computed. It is tailored for the case of solving Ax = b, x = A _1 b > 0. 

Theorem 4. (Sensitivity Theorem [17J) Let 

A = min {||p|| : p e conv ({ai,...,On})}, (21) 

p = max{||oi||,...,||a B || ) ||6||}. (22) 
Let Aq be any number such that < A < Ao- Suppose e € (0, 1) satisfies 

e < (23) 
2p 

Suppose we have computed 

p' = (aiai H h a n a n ) - a n+1 b € conv{{a 1 , . . . ,a„, ~b}) (24) 



satisfying 
Let 

Then, xq > 0, and if 
we have 

\\Ax - b\\ < e'p, (28) 
i.e. Xq is an e' -approximate solution of Ax — b. □ 
Now consider the following Two-Phase Triangle Algorithm: 



Ib'll < (25) 

XQ = ( ) . (26) 

Oin+1 «ri+l 



e ' = 2 ( 1+ Af) e ' (2?) 



Two-Phase Triangle Algorithm (A — [a\, . . . ,«„], b) 

• Phase 1. Call Triangle Algorithm({ai, . . . , a n }, 0) to get a witness p 1 E 
conv({ai, . . . , a n }). 

• Phase 2. Starting with p' in Phase 1, call Triangle Algorithm({ai, . . . , a n , —b}, 0). 



The first phase of the algorithm attempts to find a witness p' £ conv({a\, . . . , a n }) that proves is not 
in this convex hull. Any such witness p' according to ([6]) gives rise to a lower bound to Ao which in turn 
can be used to select e, see (|2"5)l in Theorem 2J From this and the complexity result in Theorem [31 we have: 

Theorem 5. Given any cq € (0, 1), in order to compute an cq- approximate solution (i.e. a solution xq > 
such that \\Axq — b\\ < eop) it suffices to set A = 0.5c?(p', 0), where p' is the witness computed in Phase 
1 of the Two-Phase Triangle Algorithm. Then in Phase 2 of the algorithm it suffices to compute a point 
p' G conv({ai, . . . ,a n , —6} so that 

Ib'll < ep, (29) 



where e satisfies 

2 —\p> (A' + \\b\\) 



£ <^ min (I, 60 y (so) 
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In particular, since Ao < p, it suffices to pick 

e < ^e . (31) 

Then the number of iterations in Phase 2 of the algorithm, K e , each of cost 0(n 2 ) arithmetic operations, 
satisfies 

Remark 2. The complexity of Phase 1 is dominated by that of Phase 2. To see this note that from Theorem 
|3]to compute a witness in testing if lies in conv({a\, . . . , a„}), the algorithm requires if a iterations. Since 
e < Aq/2/3, and So < p, it follows that 

ln§<ln-. (33) 
Thus Ka is bounded above by the bound on K e in Theorem [5l 

Remark 3. In practice, with the goal of computing an e-approximate solution, we do not need to compute 
an a priori estimate A' Q on Ao to be used in Theorem [5] We could forgo Phase 1 altogether and simply run 
Phase 2. Assuming that A is invertible, in each iteration we compute from the current iterate p' the current 
approximate solution xq, then check if it is an e -approximate solution of Ax = b. 

The following result, while not practical, suggests a lower bound to Ao can be stated based on the smallest 
eigenvalue of A T A. 



Proposition 3. Let Aq be as in (21\) . Let A m i n be the minimum eigenvalue of the matrix Q = A T A. Th 



en 



Aq = min{||AB|| : £x i = l, Xi > 0} > -^A min . (34) 



i=l 



Proof. Suppose Ao = |j4^'||, where X^Li x i = 1> x ' — ®- ^ ^ s eas Y to prove |ja;'|j > l/y/n. We thus have 

\\Ax>\\* = x' T Qx> = \\xT(^fQ(^) > ^AL„. (35) 

□ 

Before considering the general case of solving Ax — b we consider a small example. 
Example 1. Consider the 2x2 linear system. 



1 1 s - ri 



(36) 



Its solution is x — (1,2). Thus the Two-Phase Triangle Algorithm can compute an approximate solution 
to any prescribed accuracy. We will consider one iteration, skipping Phase 1. To compute an approximate 
solution we test via the Triangle Algorithm if lies in the convex hull of the set 

^=($)^={-f)^={^\ (37) 

Let the initial iterate p 1 be selected as the center of mass, p' = (a\ + — 6)/3 = (2/3, — 1/3) T . Thus 
the initial approximate solution to Ax = b is xq = {cti/ 0.2,0.2/ o/-3) T = (1, 1) T - The corresponding error is 
\\Axq — &|| = \fh. In Step 1 of the Triangle Algorithm we select 02 as the pivot point since from Remark [T] 
we have 

p' T a 2 < 0. (38) 
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The corresponding step size is 

a ~ \\a 2 -pT 4" m 

Thus 

P" = (1 - + = ^oi + -aa - i 6. (40) 

This implies the next approximate solution to Ax — b is x\ — (1,2) (using that the new coefficients are 
a\ — 1/4, Q!2 = 1/2, «3 = 1/4) which is indeed the exact solution. 

3.2 Solving A General Linear System Via The Triangle Algorithm 

In this section we consider the general case of solving Ax = b with A an invertible matrix, where it is not 
known if the solution x = A~ x b is nonnegative. Let e = (1, . . . , 1) T € K m . Then there exists t > such that 
if x is the solution to 

A(x - te) = b, (41) 

then x > 0. Thus if we let u — Ae, then 

Ax = b + tu, x>0, (42) 

is solvable. Since A is invertible ii / 0. Let 

U = min {t : A(x - te) = b, x > 0}. (43) 

If a value £ is known such that t > t*, then we can apply the Two-Phase Triangle Algorithm to solve 
(|4"2"j) . In the complexity analysis we use bounds that also depend upon t. Specifically, set 

b(t) =b + tu. (44) 

Then we may restate Theorem 2] as well as the complexity result, Theorem [SJ with ||b|| and p replaced with 
\\b(t)\\ and p(t) defined as 

/ j(t) = max{||a 1 ||,...,||a n ||,||6(t)||}. (45) 
For any given t>0we can state the following complexity in testing the solvability of (|42|) . 



Theorem 6. Given eo > 0, and any t > 0, and any lower bound A' , < Aq < Ao, the Two-Phase Triangle 
Algorithm in at most 

48pft) 2 
£ o Aq 2 

iterations, each of cost 0(n 2 ) arithmetic operations, either determines that Ax = b + tu, x > is infeasible, 
or computes Xq satisfying 

\\Ax - b\\ < e p. 



Next we show how to actually find an upper bound on f* defined in (|43j) . 

Theorem 7. Consider Ax = b and assume A is invertible. Let x* = A~ 1 b and x* be its i-th component. 
Let Q = A T A, w = A T b, and let the i-th column of Q be denoted by qi. Let g m ; n = min {||<7j||,.7 = 1, . . . ,n}. 
Let A m i n be the minimum eigenvalue of Q. Set 

Then 

x*i > > -t». (47) 

In particular, r* > > . 
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Proof. We have Qx* = w. By the Cramer's rule we have 

* det(Qi) 



1 det(Q) ' 

where Qi is the matrix that replaces column i of Q with the vector w. Thus, 



(48) 



By Haddamard's inequality on determinants we have 

\^m<(flM\)^<(flM)^- (50) 

On the other hand, det(Q) = YVj=i where Ai, . . . , A n is the set of eigenvalues of Q. Since Q is positive 
definite, these are positive and we have det(Q) > A^ in . These imply (|4"T1) . Since A((x* + r^e) — r^e) = b, 
from (|43|) it follows that + r*e > 0. Hence, > t*. But also t* > r*. □ 



3.3 An Incremental Triangle Algorithm For Solving a Linear System 

In this section we consider solving Ax = b without using an a priori upper bound t on the value of t* needed 
to make the solution of Ax = b + tu nonnegative (see (|42|l ). That is, rather than selecting a specific value for 
t > t* , we gradually increase its value (starting from £ = 0) and solve the corresponding system Ax = b + tu 
via the Triangle Algorithm to get an e-approximate solution. If for a particular value of t the Triangle 
Algorithm fails to generate an e-approximate solution we then increase t and repeat the process. However, 
in doing so we make use of the duality implied by Theorem [1] i.e. the existence of a witness computed via 
the Triangle Algorithm. 

More specifically, the Incremental Triangle Algorithm works as follows. Assume that for a given to > 
(initially set be to be zero) we have attempted to compute an e-approximate solution for Ax = b, i.e. a 
vector xq > such that 

\\A(x - t e) - 6|| < emax{||ai||, . . . , ||a„||, ||6|| }. 

If this is possible we are done. If not, then convi\\a\ 1 . . . , a n , — b(to)\) does not contain the origin, where 
b(to) = b + toil. Thus by Theorem Q] the Triangle Algorithm computes a witness, i.e. 

p'(t Q ) e conv({ai, ...,a n , -b(t )}) (51) 

such that the following set oi n+1 strict inequalities are satisfied: 

|b'(t )-ai|| < |N|, Vi = l,...,», (52) 

and 

b'(to) + &(*o)ll <||6(*o)||. (53) 
Equivalently, after expanding and simplifying (|52j) and (|53j) we get 

\\ P '(t )\\ 2 -2 P '(t Q ) T ai <0, V» = l,...,n, (54) 

\\ P '(t )\\ 2 + 2p'(t ) T b(t ) <0. (55) 

From (1511) we have 

n n+1 

P'(~to) = ^ a i a i - «n+l(ft + tpu), ^ Qij = 1, C*i > 0, Vi (56) 
t=l j=l 
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Letting 

we may write 
Thus, 

For each t define 

For i = 1, . . . , n, define 

Also, define 



p — 22 a i a i - a-n+ib, 
»=i 

p'(t ) =p' - t Q a n+ iu. 
p'(t ) T a t = p' T a t - t a n+ iu T ai. 
p'(t) =p' - ta n+1 u. 
g l (t) = \\p'(t)\\ 2 -2p'(t) T a l . 
g n+1 (t) = \\p'(t)\\ 2 + 2p'(t) T b(t). 



(57) 

(58) 
(59) 
(60) 
(61) 
(62) 



It is easy to verify that for i = 1, . . . , n we have 

9i(t) = t 2 a 2 n+1 \\u\\ 2 - 2ta n+1 {p' - ai ) T u + \\p'\\ 2 - 2p a \. (63) 
The coefficient of t 2 in g n +i(t) can be shown to be 

a„+iK+i -2)||u|| 2 . (64) 



We now prove 

Proposition 4. Suppose a n +i ^ 0. For each i — l,...,n, we can select a real number t{ > to such that 
9i(ti) > 0. Moreover, g n +i(t) < for all t. 

Proof. Consider a quadratic polynomial q(t) = C2t 2 + c\t + co, where C2 > 0. We claim if q(t) takes a negative 
value at some point to, then it must have a real root t' > to- The proof of this follows from the fact that 
as t approaches infinity, q(t) approaches infinity so that q(t) changes sign from negative at to to a positive 
value. Thus q(t) must be zero at some t' > to. For each % = 1, . . . , n, the coefficient of the quadratic term in 
gi(t) is positive. This follows from the fact that u ^ and a n +i ^ 0, see (|63]). Hence applying the argument 
on q(t) we conclude that each 9i(t), i = 1, . . . , n, has a root ti larger than to- 

Since the coefficient of the quadratic term in g n+ i(t) is a n +i(a„+i — 2)||u|| 2 and < a n +i < 1, the 
coefficient is negative. Since g n+ i(to) < 0, it follows that g n +i{t) < for all t. □ 

Let t' = min{ti, . . . ,t n }. Increasing to to t' , at least for one i — 1, . . . ,n we will have gi(t' ) > 0. Once t' 
is computed we have a pivot point and we reapply the Triangle Algorithm, testing if Ax — b + t' u, x > is 
feasible. In doing so we can start the Triangle Algorithm with p'(t' ) = p' + t' b. In practice we do not need 
to find the exact value of ti which requires solving the quadratic equations <fe(t) = 0, i = 1, . . . , n. We can 
simply compute an upper bound on the roots. Such upper bounds can be computed easily. Additionally, we 
can choose t' to be such that t' — to differs by a natural number, guaranteeing that this iterative process 
would eventually select a value t' exceeding the value t* , see (1431) . A formal description of the Incremental 
Triangle Algorithm is as follows. 
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Incremental Triangle Algorithm^, b) 

• Step 0. (Initialization) Let u = Ae, e = (1, . . . , 1) T e K m . Select p' — Aa - a n+1 b e 
conv(S), where a = (oti, . . . , a n ) T £ R", J27=i a i = 1> a « — 0- Set to = 0. 

• Step 1. Given p' = Aa — a n+ \b € conv(S), set xo = a/a n +i. Define To according to 

£J(t ) = ||Ax - (6 + r w) || = mini \\Ax - (b + tu)\\ : * > *o V- (65) 



Replace to with to. If £(to) < set x' = (xq — toe), stop. 

• Step 2. If p'ito) = p' — a n+ itou is a witness with respect to conv({ai, . . . , a„, — (fe + to u )}): 
go to Step 3. Otherwise, call Triangle Algorithm ({ax, ■ ■ ■ , a n , — (b + t u)}, 0) to compute a 
new iterate p"(to) (see (fTU|) in Triangle Algorithm): 

p"(*o) =p" - Pn+itou, where p" = A/3 - /3 n+ ib € conv({a\, ...,«„, -(& + to«)}))- (66) 

Replace p' with p", a with /3, and a n +i with /3„+i. Go to Step 1. 

• Step 3. Compute t' , the smallest value t such that <&(<) > for some i = 1, . . . , n (see 
(J63l)). Replace t with t . Go to Step 2. 



Remark 4. The computation of To in Step 1 is an auxiliary step to improve the error for a given xq. For a 
given e, eventually x' will give the desired e-approximate solution. 

We will consider an example. In the example below we will implement the algorithm with and without 
the auxiliary Step 1 to show different scenarios. 

Example 2. Consider the 2x2 linear system (I57|) . having x = (—1, — 2) T as its solution. Thus f* = 2. We 
consider one iteration in solving Ax = b by the Incremental Triangle Algorithm. 



J ft) " I- 

We have it = Ae — (1, 2) T . We set to = and test if is in the convex hull of the set 



(67) 



v fli= (9'* = ("i)'- 6= (3 )}- (68) 

Suppose we select p' = (0,3/2) T . It is easy to check that 

p' = a\a\ + a2<ii — a^b, a\ = 1/4, o?2 = 1/2, 0:3 = 1/4. (69) 
Thus the initial approximate solution to Ax = b is 

*o = (^) r = (l,2f. (70) 

a 3 a 3 

To compute To we minimize 

B(t) = ||Auo-(6 + t«)|| = ||(-t,6-2t) T || (71) 

for t > to. Note that E(to) — 6. However, the minimum of E(t) is attained at To = 12/5 = 2.4 and 
E{t ) — 6/V5 ~ 2.68. For simplicity of the computation we round the value of t to 2. Next, we replace 
t with tq. We have, E(2) = 2v / 2- In Step 2 of the algorithm we have p'(to) = p' — a 3 t u = (—1/2, 1/2) T . 
From Remark [T] we may select a\ as the pivot point as we have, 

p'(t ) T ai = -1/2 < 0. (72) 
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(73) 



The corresponding step size is 

p'{t Q ) T ( ai -p'{t )) _ 2 
\\ ai - P '(t W 13" 
Thus 

,,, x / 2 , , 2 11,1 1 1„ 2 

/ '(to) = (1 - + Y3 fll = 13 ^I ai + 2° 2 ~ 4 (6 + n)) + 13 ai = 

19 11 11,, „ x 

52 ai + 26 a2 -52 (6 + 2M) - (74) 
This replaces p'(to), and the new p' becomes 

P = 52 ai + 26* - 52^ (75) 

The next approximation to the linear system is xq — (19/11, 2) T . The corresponding error, E(to) = 
V936/11 ~ 2.78 which is less than 2y2- Now we optimize again by computing tq as the minimum of 
E(t) = \\Axq — (b + tu)\\. The minimum occurs for tq = 164/55 and is approximately 1.7. We replace to 
with this and repeat Step 2. 

Rather than continuing with this iteration, we will next consider the same example but without computing 
To- This will allow us to implement one iteration of Step 3 in this example. Starting again with the point 
p' = (0, 3/2) we can easily check that it is a witness, i.e. 

P ax > , p a 2 > — — , -p b> —^—- ( 76 ) 
p ' {t )=p'-ta,u=(- t - Z -- t -) T . (77) 



Then 



Using this to compute p'(t) Oj, for i = 1, 2 we get 



9l (t) = b'(t)|| 2 - 2p'{t) T ai = ^ + \t-\, (78) 

g 2 (t) = Hp'WII 2 - 2p'{t) T a 2 = ^t 2 - t - \. (79) 
Also, we have b(t) =b + tu = (0, -3) T + t(l, 2) T = (t, -3 + 2t) T . So we have 

93(f) = b'WII 2 + 2p'(tfb(t) = + Y*-i- ( 8 °) 

As analyzed for the general case, the quadratic term of ga(t) has a negative coefficient so that ga(t) remains 
negative for all t. We select t' to be the positive root of gi{t), namely (— 8 + v / 304)/10 « .94. The positive 
root of g 2 (t) is larger than this value. 

The algorithm then replaces to with t' , and moves to Step 2, setting p'(io) = p' + toot^u, then checking 
if it is a witness with respect to conv{a\ 1 a 2 , — b(to)). Note that since = 2, we expect that if p'(to) is not 
a witness, subsequent iterations of the algorithm (without the optimization step that computes to) must 
eventually produce a witness. 

Remark 5. The Incremental Triangle Algorithm avoids using an a priori estimate for £*. It increases t 
incrementally, starting with t — 0, generating a sequence of approximate solutions Xk that converge to the 
solution of Ax = b. Here Xk and Xk+i may correspond to the same value of t, or to two consecutive values 
of t. The increase in the t value in the incremental algorithm is done in a conservative fashion. In order 
to guarantee that it eventually gets to be large enough, e.g. at least as large as (see (|43]) ). or large 
enough such that an e-approximate solution is possible, we will make sure that the difference between two 
consecutive t values is at least a natural number no- 
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Another possibility in increasing the t value after each report of the infeasibility of Ax — b(t), x > 
via a witness, is to double it and add one. Clearly, this would require at most 0(log 2 (t* + 1)) calls to the 
Triangle Algorithm. However, the potential advantage in the Incremental Triangle Algorithm is that once for 
a particular value t we get a witness p' , for the next value, say t' , if still not sufficiently large, the algorithm 
will find a new witness p" quickly, in a few iterations. 

4 Final Remarks 

In this paper we have proposed two novel iterative methods for solving linear system of equations approxi- 
mately to within prescribed errors. These approximation algorithms are based on the Triangle Algorithm for 
a convex hull problem, |17] . Undoubtedly testing the practical performance of these algorithms will require 
extensive experimentation, as well as comparisons with existing computational results via exact or iterative 
algorithms. Ideally, these experimentations should be applied to different types of matrices, from large to 
sparse matrices. We hope to carry out some experimentation to assess the performance of the Triangle Al- 
gorithm for the convex hull problem, as well as for solving a linear system of equations via the two proposed 
algorithms. Optimistically, the Triangle Algorithm will perform well in computing approximate solutions to 
linear systems, perhaps with much better practical performance than its theoretical worst-case complexity. 
Potentially, the proposed algorithms can also be combined with the existing iterative algorithms for linear 
systems. The simplicity of the Triangle Algorithm and its theoretical performance raises the optimism that 
it will encourage new research and applications. 
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